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A new numerical finite difference code has been devel- 
oped to solve a thermal convection of a Boussinesq fluid 
with infinite Prandtl number in a three-dimensional spheri- 
cal shell. A kind of the overset (Chimera) grid named "Yin- 
Yang grid" is used for the spatial discretization. The grid 
naturally avoids the pole problems which are inevitable in 
the latitude-longitude grids. The code is applied to numeri- 
cal simulations of mantle convection with uniform and vari- 
able viscosity. The validity of the Yin- Yang grid for the 
mantle convection simulation is confirmed. 



1. Introduction 

From the middle of 1980s, numerical simulation codes 
for the thermal convection with infinite Prandtl number in 
three-dimensional (3-D) spherical shells have been developed 
to solve the mantle convection of terrestrial planets. The dis- 
cretization methods employed in these codes can be divided 
into three categories; the spectral method [Machetel et al, 
1986; Glatzmaier, 1988; Bercovici et al, 1989; Zhang and 
Yuen, 1995; Harder and Christensen, 1996], the finite ele- 
ment (FE) method [Baumgardner, 1985; Bunge and Baum- 
gardner, 1995; Zhong et al., 2000; Tabata and Suzuki, 2000; 
Richards et al., 2001], and the finite volume (FV) method 
[Ratcliff et al., 1996; Iwase, 1996]. The spectral method, 
which can be an effective method for spherical fiows [e.g., 
Fornberg, 1996; Fornberg and Merrill, 1997], had found to 
be unsuitable to mantle convection simulations because of 
intense spatial variation of the viscosity of mantle rock. A 
new method based on multilevel wavelet algorithm [Vasi- 
lyev et al, 1997] can treat the spatially localized physical 
properties and has a great potential usefulness in mantle 
convection simulations. Its application to a spherical shell 
model is, however, still remains a challenging task. Among 
the grid-based FE, FV and finite difference (FD) schemes, 
the FV and FD methods are more desirable than FE for mas- 
sively parallel vector computers because of their feasibility 
of optimization. Another advantage of the FD method is 
its fiexibility; the extension to higher-order schemes, which 
might be important to obtain accurate solutions of thermal 
convection with very large Rayleigh numbers [e.g., Larsen 
et al., 1997], is relatively easy. 

One of the most popular computational grids in the spher- 
ical polar coordinates (r, 6, <p) is latitude-longitude {6, (j))- 
grid, which is defined by intersections of latitude and lon- 
gitude circles on a sphere (Fig. la). It is widely recog- 
nized that the {6, 4))-gnd has the "pole problems" that 
refer to two different kinds of difficulty in numerical cal- 
culations; one is the coordinate singularity on the poles 
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(6^ = 0,7r); and the other is the grid convergence near the 
poles. The pole problems have been considered as serious 
difficulties in the community of mantle convection simula- 
tion. To avoid the coordinate singularity, special cares have 
to been taken. In the FV method, for example, all the 
physical variables are arranged not to reside on the pole 
grids [Ratcliff et al., 1996; Iwase, 1996]. The problems 
of the grid convergence is more serious than the coordi- 
nate singularity: It causes not only the grid redundancy, 
but also the severe restriction on the time-step due to the 
Courant-Friendrichs-Levy (CFL) condition. In the {9, (p)- 
coordinates, the grid spacing on the spherical surfaces is 
extremely non-uniform as Fig. la shows. The largest grid 
spacing AX is given in the equator; AX = 2tt/N^, where 
A'^^ is the grid number in the (^-direction, while the small- 
est grid spacing Aa; is given at the nearest latitude to the 
poles; Aa; = rsin(7r/Xe) x (2n/N^) ~ 2n^r/{NeN^), where 
Ng is the grid number in the ^-direction. So the ratio 
AX/ Ax ~ Ne/iT increases in proportional to the grid imm- 
ber. This means that the time-step restriction becomes ex- 
tremely severe for large scale simulations with fine grids. To 
avoid the impractically small time-step, one has to invoke 
quasi-uniform grid spacing over the sphere. The FE based 
codes referred above employed carefully designed grid cells 
for that purpose. For example, a FE mantle convection code 
named CitcomS has nearly uniform resolution in both polar 
and equatorial regions [Zhong et al., 2000]. However, a FD 
or FV based mantle convection code that overcomes both of 
the pole singularity and the grid convergence have not been 
reported so far. 

Here we employ a new grid system for spherical shell ge- 
ometry, named "Yin- Yang grid", which has been proposed 
recently by Kageyama and Sato [2004]. The Yin- Yang grid 
is composed of two component grids that have exactly the 
same shape and size (Fig. lb). They partially overlap each 
other on their boundaries (Fig. Ic). Following the over- 
set (Chimera) grid method [Chcsshire and Hcnshaw, 1990], 
data on the boundaries of the component grids are matched 
by interpolation. A component grid of the Yin- Yang grid 
is actually a low latitude part of the {6, </))-grid. As it is 
apparent in Fig. lb, the Yin- Yang grid has neither a coor- 
dinate singularity, nor grid convergence; the grid spacings 
are quasi-uniforms on the sphere (see Kageyama and Sato 
[2004] for more details on this grid). 

In this paper, we apply the Yin- Yang grid for the numeri- 
cal simulation of mantle convection. To confirm the validity 
of the Yin- Yang grid, we have performed benchmark tests 
with published numerical codes for steady convections. We 
also apply the Ying-Yang grid for time-dependent mantle 
convections with uniform and variable viscosity. 

2. Model and Numerical IVIethods 

We model the mantle convection as a thermal convection 
of a Boussinesq fiuid with infinite Prandtl number heated 
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from bottom of a spherical shell. The ratio of the inner ra- 
dius (r = ro) and the outer radius (r = ri) is 0.55. The 
normalization factors for the non-dimensionalization of the 
length, velocity, time and temperature are f i = 6371 km (the 
Earth's radius), A/ri, ri/k and AT = T^ot — Ttop, respec- 
tively, where ft is the thermal difFusivity, and Ttot and Ttop 
are the temperatures on the bottom and top surfaces. The 
hat stands for dimensional quantity. The non-dimensional 
equations of mass, momentum, and energy conservation gov- 
erning the thermal convection are. 



V • V = 0, 

-Vp -I- V ■ {r]e) + RaCTer, 
dtT = V^T - V ■ VT, 



(2) 
(3) 



where v is the velocity vector, p the dynamic pressure, 
r the temperature, t the time, e the strain-rate tensor, 
and Br is the unit vector in the r-direction. The con- 
stant parameter ^ is (d/fi)'^ = 0.45'\ where d is the 
thickness of the shell, 2890 km (the Earth's mantle). We 
assume that viscosity r? depends only on temperature; 
rj{T) = rjref &x-p[—E {T — Tref)], where T^e/ is the refer- 
ence temperature, and rjref is the reference viscosity at 
Tref- The parameter E denotes the degree of viscosity con- 
trast between the top and bottom surfaces. The viscos- 
ity contrast across the spherical shell is defined by 7,, = 
ri{Ttop)/ri{Tbot) ~ exp(i5). The Rayleigh number is defined 
by Ra = pgaATcr / kfiref , where p is the density, g the 
gravitational acceleration, and a is the thermal expansivity. 
The mechanical boundary conditions at the top and bottom 
surface are immpermiable and stress-free. The boundary 
conditions for T are Tbot = 1 and Ttop = 0. 

Wo use the collocated grid method [e.g., Ferziger and 
Peric, 2002]; all the primitive variables, v, p and T, are de- 
fined on the same grid points. Equations (l)-(3) arc solved 
by the FD discretization with second-order accuracy. The 
SIMPLER algorithm [Patankar, 1980; Ferziger and Peric, 
2002] is applied to solve v and p from eqs. (1) and (2). 
The Crank-Nicolson method is used in eq. (3) for the time 
stepping. The upwind difference method is applied for the 
advection term in eq. (3). With the Yin- Yang grid method, 
we simultaneously solve eqs. (l)-(3) for each component grid. 
We use a successive over-relaxation (SOR) method as the it- 
erative solver required in the SIMPLER algorithm and the 
energy equation. The horizontal boundary values of each 
component grid are determined by linear interpolation from 
the other component grid. The interpolation is taken at 
each SOR iteration. (We confirmed that the interpolation 
procedure has no numerical mischief on the calculations.) 
The grid size is 102 x 102 x 204 (in r-, 6-, and </>-directions) . 
Wo have confirmed that this size is enough to resolve all 
the convections studied in this paper. Time development of 
the convection is calculated until averaged quantities, such 
as Nusselt number and root-mean-square velocity, become 
stationary. 

3. Benchmark Tests 

The thermal convection in the spherical shell with in- 
finite Praudtl number has two stable solutions with poly- 
hedral symmetry when the Rayleigh number is low [e.g., 
Schubert et al., 2001]. The two solutions are found by lin- 
ear theory [Busse, 1975; Busse and Riahi, 1982] and con- 
firmed by numerical simulations [Bercovici et al., 1989; Rat- 
cliff et al, 1996]: One solution is a convection with the 
tetrahedral symmetry which has four upwellings; the other 



has the cubic symmetry with six upwellings. To confirm 
these symmetric solutions and their stabilities, we performed 
two simulations with different initial conditions of tem- 
perature field; T{r,e,(j)) = T^ondir) + Tprtbir,9,(l>), where 
Tcondir) = ro(ri — r)/r(ri — ro) is the purely conductive 
profile, i.e., V Tcond(r) = 0, with the thermal boundary 
conditions given above. The perturbation term Tprtb{r, 9, 4>) 
is given by. 



Tprtb{r,0,<j)) = A \Yi{e,,j)) + n{e,,p)\ simr{r-ro), 



(4) 



^j^'j for the tetrahedral symmetric solution, and 



Tprtbir, e,<j>) = x \Y4°{e, 4>) + ^Y^\e, <t>) + n{6, ^ 

X sin7r(r — ro), (5) 

for the cubic symmetric solution, where Yi"{0, (p) is the fully 
normalized spherical harmonic functions of degree i and or- 
der rn. The Yi"^{6, (f>) terms in eqs. (4) and (5) determine the 
solution's symmetry. The other term 0.{6, 0) is for secondary 

perturbation. We set ^{6, <t>) = oj y!^^o ^'!"(^' 

We have performed benchmark tests with published nu- 
merical mantle convection codes that employed various nu- 
merical schemes. Following Richards et al. [2001] and 
Ratcliff et al. [1996], we performed simulations of uni- 
form (7^ = 1) and variable (7,, = 20) viscosity convec- 
tions with both the tetrahedral and cubic steady symme- 
tries when A — 10~^ and a; = (i.e., no secondary per- 
turbations). The Rayleigh number Rai/2 is defined by the 
reference viscosity rfref at Trcj = 0.5 [Ratcliff et al., 1996]. 
Nusselt number at the surface and root-mcan-square ve- 
locity of entire domain were calculated on convections at 
Rai/2 = 2.0 X 10^ ~ 1.4 X 10". The resuhs of the benchmark 
tests are summarized in Table 1. In spite of the differences of 
the discretization methods, numerical techniques, and num- 
ber of grid points among the codes, we found that the results 
from our code agree well with them within a few percent or 
even better and confirmed the validity of our code. 

4. Unsteady Convection Problems 

The steady convections become time-dependent when the 
Rayleigh number is increased. Since the Earth's mantle 
is obviously time-dependent convection with high Rayleigh 
rmmbcr, the transition of convection from steady to un- 
steady state is important. We tried a series of simulations 
with various Rabot (the Rayleigh number defined by the ref- 
erence viscosity at the bottom surface, i.e., Ti-e/ = Tbot) from 
the critical number for convection onset (« 712) [Ratcliff et 
al., 1996] to 10'. The perturbation amplitudes A and ui are 
taken to be lO"'^ and 10~^, respectively. Shown in Fig. 2 
are the iso-surfaces of temperature at Rubot = 10* and 10^ 
after 200,000 time-steps. Figure 2a and 2b indicate that, at 
Rabot ~ iC, the convection patterns are in steady states, 
maintaining each symmetry, in spite of the existence of the 
secondary perturbations in the initial conditions. This is 
consistent with earlier results [Bercovici et al., 1989; Rat- 
cliff et al., 1996] in which the secondary perturbation was 
not explicitly imposed, i.e., to = 0, though. 

When Rabot ~ 10"", the convection patterns become 
weakly time-dependent. The geometrical symmetry in this 
Rayleigh number is broken. This disagrees with the result of 
Ratcliff et al. [1996]. Notice that, in the right panel of Fig. 
2b, all the six upwelling plumes have the same diameters in 
our results. The corresponding case by i?atefejff ei ai. [1996], 
in which a FV scheme on the (O,0)-grid is used, shows a 
symmetric pattern about equator and appears to remain in 



a steady state [of. RatchjJ et ai, 1996, Fig. 6]. These obser- 
vations suggest that the low Rayleigh number convections 
around Ratot ~ 10'' are numerically affected by coordinate 
singularity and the grid convergence in the (6, (f>)-gnd. On 
the other hand, the pole effects are removed in our code by 
making use ol the Yin- Yang grid. 

It is known that variable viscosity with strong tempera- 
ture dependence induces drastic effects on the convection 
pattern in 3-D Cartesian model with large aspect ratio 
and also in the spherical shell model [RatchjJ et al., 1997; 
Trompert and Hansen, 1998]. To confirm this effect in our 
model, we performed simulations with variable viscosity. 
Taking eq. (4) as the initial temperature perturbation, we 
first calculated an isoviscous convection at Rabat = 10®. The 
obtained solution, which is shown in Fig. 3a, is strongly 
time-dependent and exhibits complex feature in contrast to 
the case at Ruhot = 10^ (the right panel of Fig. 2a) . We grad- 
ually increased 7^ from 1 (isoviscous case) up to 10^. We 
obtained a convection regime that has cold and rather thick 
thermal boundary layer on the top surfaces (Fig. 3b). The 
large aspect ratio of convecting cells in this regime is consis- 
tent with the previous results obtained by the 3-D Cartesian 
model with large aspect ratio as well as spherical shell model 
with moderately strong temperature-dependence of viscos- 
ity (7r) ~ 10^) [Ratcliff et ai, 1997]. Our results show that 
the underlying convection patterns with larger aspect ratio 
of degrec-2 come to dominate. The two cells structure that 
consists of one sheet-like downwelling along a great circle of 
spherical shell and two mushroom-shaped upwelling plumes 
is formed. 

5. Conclusions and Discussion 

We have developed a new numerical simulation code to 
solve the thermal convection of a Boussinesq fluid with in- 
finite Prandtl number using a second-order FD method on 
newly devised spherical overset grid named Yin- Yang grid. 
The validity of the Yin- Yang grid for the mantle convec- 
tion simulation is confirmed by benchmark tests. Our code 
is powerful and unique FD based code that can solve both 
the uniform and the strongly variable viscosity convections. 
The Yin- Yang grid is suitable to solve the mantle convection 
problems because it automatically avoids the pole problems 
that arc inevitable on the {6, (j))-grid. In the isoviscous case 
with cubic symmetry at Ratot = 10'', the convection pattern 
has a weak time-dcpondoncc in our Yin- Yang grid, while it 
was steady with strange asymmetry of the plume sizes be- 
tween those on the poles and those in the equator in the pre- 
vious FV scheme on the {6, (?i)-grid. This discrepancy might 
be a consequence of the grid convergence near poles in the 
{9, 0)-grid. Our result implies that large-scale (low degree) 
convcctive structures are easily affected numerically by the 
poles when {9, (f>)-grid is employed. The quadrulpolc convec- 
tion patterns is obtained when large viscosity contrast with 
three orders of magnitude is introduced when Rabat = 10®. 

To follow mantle convection for geophysical time-scale 
(~10** years), the computational time-step A.t is critically 
important in rmmcrical sinmlations. As we described in sec- 
tion 1, the time-step is determined by the CFL condition 
by the smallest grid spacing. For {9, 0)-grid, Aa;(= Aa;e^) 
is determined by the azimuthal grid spacing at the nearest 
grids to the pole. On the other hand for the Yin- Yang grid, 
Aa;(= A.xyy) is determined by the azimuthal grid spacing 
at 9 = -k/A (or 37r/4). Therefore the ratio of time-steps Vt 
between two grids is. 

Ft oc Aa;e<^/Aa;i'y = sin(7r/Are)/sin(7r/4) w lAn/Ne. (6) 

Taking Ne = 102 as employed in this paper. Ft « 0.04. This 
means that the total computational time is significantly re- 
duced by the factor of 1/25 by making use of the Yin- Yang 
grid. 
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Figure 1. The latitude-longitude {9,(j))-grid and new 
spherical overset grid named "Yin- Yang grid", (a) The 
(8, </))-grid. (b) Two component grids of the Yin- Yang 
grid. They are identical (same shape and size); the low 
latitude part {n/4 < 6 < 3n/4, -Sn/A < 4> < 37r/4) of 
the (9, 0)-grid. (c) They partially overlap each other at 
their interface to cover a spherical surface in pair (see 
text and Kageyama and Sato [2004] for details). 



(a) Tetrahedral 

Rabat =10'' Rabot=10^ 




Figure 2. The iso-surface renderings of temperature 
(T = 0.4) started from the initial conditions of (a) tetra- 
hedral, and (b) cubic symmetries. The left and right 
panels on each figure show the cases at Ratot ~ 10^ and 
10^, respectively. 



(a) , (b) 



Figure 3. The iso-surfacc renderings of residual temper- 
ature ST (i.e., the deviation from horizontaUy averaged 
temperature at each depth) for the cases at 7^ = (a) lO", 
and (b) 10^. Blue iso-surfaces stand for 6T of (a) —0.10, 
and (b) -0.25. Yellow iso-surfaces for 6T of (a) -1-0.10, 
and (b) -f 0.25. Red spheres indicate the bottom surface 
of spherical shell. 



Table 1. The benchmark test of Nusselt numbers at the top surface {Nu) and RMS velocity (Vrms) of the entire mantle" 
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a ur^jQ,, (jQjjQi^Qg tctrahcdral ("T") or cubic ("C") symmetric solutions. The abbreviated code names "G188" is for Glatzmaier 
[1988], "Br89" Bercovici et al. [1989], "HC96" Harder and Christensen [1996], "Rt96" RatclijJ et al. [1996], "Iw96" Iwase [1996], 
"ZhOO" Zhong et al. [2000], "TSOO" Tabata and Suzuki [2000], "ReOl" Richards et al. [2001], and "YK04" is for our code. The "SP" 
in parentheses under each code name denotes spectral method, and see text for "FV", "FE" and "FD". (Note that, in this benchmark 
test, the normalization factor used to non-dimensionalize the length is the Earth's radius fi, not d.) 



